Actividad 1 M4

1. Construye un modelo de dinámica de sistemas basado en el caso anterior. Calcula el valor de los exponenetes y las variables de flujo omitidas de tal manera que el colapso de la población Maya ocurra en el año 800. Simula 2000 años de evolución iniciando en el año 1000 A.C

DIAGRAMA CASUAL DEL CASO

library("deSolve")

#La variable exógena growth population rate se calculó utlizando la siguiente fórmula: 
#Anual population rate= 1-(final value/starting value)^1/n
#Donde: 
#final value = 2 
#starting value = 1
#n = 408 años
#1-(2^(1/408))-1 = 0.0017
parameters<-c(
growth.population.rate=0.0017, #[% per year] 
emigration.ratio=0.05, #[% per year]
consumed.food.per.person=400, #[kg per person per year]
intensity.of.agriculture=1, 
x=1.3595 #Valor calculado, ver NOTA 1. 
)

InitialConditions <- c(
  population=100000, #[persons]
  fertility.of.agricultural.land=5000000, #[kg/km2]
  agricultural.land=8, #[km2]
  forest.land=5000 #[km2]
  )

#-1000: El signo negativo representa el tiempo B.C, por lo cual es el año 1000 B.C.
#2000: El signo posistivo representa el tiempo A.D, por lo cual es el año 2000 A.D.
times <- seq(-1000,2000,1)

#Definimos método de integración
 intg.method<-c("rk4")

#Especificamos modelo
MAYAS <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {

#auxiliary variables
food.produced <- agricultural.land*fertility.of.agricultural.land #[kg]
demand.food <- population*consumed.food.per.person #[kg]
food.gap <- demand.food - food.produced #[kg]
people.without.food <- food.gap/consumed.food.per.person #[persons]

#flow variables
natural.population.increase <- growth.population.rate * population #[persons]
emigration <- emigration.ratio*people.without.food #[persons]
fertility.losses <- fertility.of.agricultural.land* min(2,(agricultural.land/forest.land)^x)/intensity.of.agriculture #[kg/km2]
deforestation <-min(food.gap/max(fertility.of.agricultural.land,1),forest.land/4)/intensity.of.agriculture #[km2]
  
#state variables
dpopulation <- natural.population.increase - emigration #[persons]
dfertility.of.agricultural.land <- (-1)*fertility.losses #[kg/km2]
dagricultural.land <- deforestation #[km2]
dforest.land <- (-1)*deforestation #[km2]

list(c(dpopulation, dfertility.of.agricultural.land, dagricultural.land, dforest.land))
    })
}

out.normal <- ode(y = InitialConditions,
           times = times,
           func = MAYAS,
           parms = parameters,
           method =intg.method)

datos <- data.frame(out.normal)

NOTA 1 ¿Cómo se estimo la variable x?

2. ¿Qué pasa con la población? Gráfica el comportamiento de esta variable, ¿corresponde el comportamiento observada? En el año 800 la fertitly of agricultural land llega a cero, la población tienede a un valor extremadamente cercano a cero, esto debido a que la oferta de alimento decrece mientras que la demanda del mismo aumenta, ocasionado que crezca la población maya sin acceso a alimentación lo cual se traduce en emigración.

library(ggplot2)
#Gráfica de población 
 population.collapse <- ggplot (datos, aes(x = time, y = population)) + 
   geom_line(color="skyblue3")+
  xlab('Time (years)') + 
  ylab('Mayan Population (persons)') +
  ggtitle('Collapse of Mayan population') + 
  theme_minimal()+
  theme(plot.title = element_text(hjust = 0.5))

population.collapse

3. Introduce cambios en el modelo que hagan más realista el comportamiento. Dado el compartamiento observado en las variables de estado, el modelo presenta un comportamiendo de acuerdo a lo planteado en el caso.

#Gráfica de fertility.of.agricultural.land 
 ggplot (datos, aes(x = time, y = fertility.of.agricultural.land)) + 
   geom_line(color="darkgoldenrod2")+
  xlab('Time (years)') + 
  ylab('Fertility of agricultural land (kg/km2)') +
  ggtitle('Behaviour: Fertility of Agricultural Land') + 
  theme_minimal()+
  theme(plot.title = element_text(hjust = 0.5))

#Gráfica de agricultural.land 
 ggplot (datos, aes(x = time, y = agricultural.land)) + 
   geom_line(color="chocolate1")+
  xlab('Time (years)') + 
  ylab('Agricultural land (kg/km2*year)') +
  ggtitle('Behaviour: Agricultural Land') + 
  theme_minimal()+
  theme(plot.title = element_text(hjust = 0.5))

 #Gráfica de forest
ggplot (datos, aes(x = time, y = forest.land)) + 
   geom_line(color="darkgreen")+
  xlab('Time (years)') + 
  ylab('Forest Land (km2)') +
  ggtitle('Behaviour: Forest Land') + 
  theme_minimal()+
  theme(plot.title = element_text(hjust = 0.5))

4. ¿Qué tan sensible es el modelo a cambios marginales en los parámetros o a cambios estructurales? Cambio en growth population rate

library("deSolve")
library(ggplot2)

#Normal
parameters1<-c(
growth.population.rate=0.0017,   
emigration.ratio=0.05, 
consumed.food.per.person=400, 
intensity.of.agriculture=1,
x=1.3595
)

#Cambio marginal 1
parameters2<-c(
growth.population.rate=0.0017,  
emigration.ratio=0.08, 
consumed.food.per.person=400, 
intensity.of.agriculture=1, 
x=1.3595)

#Cambio marginal  3
parameters3<-c(
growth.population.rate=0.0017,  
emigration.ratio=0.1, 
consumed.food.per.person=400, 
intensity.of.agriculture=1,
x=1.3595)

InitialConditions <- c(
  population=100000, 
  fertility.of.agricultural.land=5000000,
  agricultural.land=8, 
  forest.land=5000 )

times <- seq(-1000,2000,1)

# Definimos método de integración
 intg.method<-c("rk4")

#Especificamos modelo
MAYAS <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {

#auxiliary variables
food.produced <- agricultural.land*fertility.of.agricultural.land
demand.food <- population*consumed.food.per.person
food.gap <- demand.food - food.produced 
people.without.food <- food.gap/consumed.food.per.person

#flow variables
natural.population.increase <- growth.population.rate * population
emigration <- emigration.ratio*people.without.food
fertility.losses <- fertility.of.agricultural.land* min(2,(agricultural.land/forest.land)^x)/intensity.of.agriculture
deforestation <-min(food.gap/max(fertility.of.agricultural.land,1),forest.land/4)/intensity.of.agriculture
  
#state variables
dpopulation <- natural.population.increase - emigration
dfertility.of.agricultural.land <- (-1)*fertility.losses
dagricultural.land <- deforestation
dforest.land <- (-1)*deforestation

list(c(dpopulation, dfertility.of.agricultural.land, dagricultural.land, dforest.land),emigration.ratio=emigration.ratio)
    })
}

out1 <- ode(y = InitialConditions,
           times = times,
           func = MAYAS,
           parms = parameters1,
           method =intg.method)

out2 <- ode(y = InitialConditions,
           times = times,
           func = MAYAS,
           parms = parameters2,
           method =intg.method)

out3 <- ode(y = InitialConditions,
           times = times,
           func = MAYAS,
           parms = parameters3,
           method =intg.method)

datos1 <- data.frame(out1)
datos2 <- data.frame(out2)
datos3 <- data.frame(out3)

#Gráfica de población 
 ggplot() +               
  geom_line(data = datos1, aes(time,population), 
             color = "darkgreen")+
   
  geom_line(data = datos2, aes(time,population), 
             color = "orange")+
   
    geom_line(data = datos3, aes(time,population), 
           color = "blue")+ 
   
  labs(x = "Time (years)", y = "Mayan population (persons)")+
  ggtitle("Mayan population") +  theme_minimal()+
  theme(plot.title = element_text(hjust = 0.5))

 #Gráfica de fertiliy of agricultural land 
 ggplot() +               
  geom_line(data = datos1, aes(time,fertility.of.agricultural.land), 
             color = "darkgreen")+
   
  geom_line(data = datos2, aes(time,fertility.of.agricultural.land), 
             color = "orange")+
   
    geom_line(data = datos3, aes(time,fertility.of.agricultural.land), 
           color = "blue")+ 
   
  labs(x = "Time (years)", y = "Fertility of agricultural land (kg/km2)")+
  ggtitle("Behaviour of fertility of agricultural land (kg/km2)") +  theme_minimal()+
  theme(plot.title = element_text(hjust = 0.5))

 #Gráfica de  agricultural land 
 ggplot() +               
  geom_line(data = datos1, aes(time,agricultural.land), 
             color = "darkgreen")+
   
  geom_line(data = datos2, aes(time,agricultural.land), 
             color = "orange")+
   
    geom_line(data = datos3, aes(time,agricultural.land), 
           color = "blue")+ 
   
  labs(x = "Time (years)", y = " Agricultural land (km2)")+
  ggtitle("Behaviour of  agricultural land (km2)") +  theme_minimal()+
  theme(plot.title = element_text(hjust = 0.5))

 #Gráfica de  forest land 
 ggplot() +               
  geom_line(data = datos1, aes(time,forest.land), 
             color = "darkgreen")+
   
  geom_line(data = datos2, aes(time,forest.land), 
             color = "orange")+
   
    geom_line(data = datos3, aes(time,forest.land), 
           color = "blue")+ 
   
  labs(x = "Time (years)", y = " Forest land (km2)")+
  ggtitle("Behaviour of forest land (km2)") +  theme_minimal()+
  theme(plot.title = element_text(hjust = 0.5))

Descripción de la gráfica: El color verde de la gráfica representa el comportamiento de las variables de estado del sistema con una growth population rate del 5%, el color naranja representa el comportamiento con una growth population rate del 8%, mientras que la línea de color azul representa el comportamiento con un 10% de growth population rate. Con base a las gráficas podemos decir que el sistema no es sensible a cambios marginales en los párametros.

Cambio marginal en el parámetro consumed.food.per.person

library("deSolve")

#Normal
parameters1<-c(
growth.population.rate=0.0017,  
emigration.ratio=0.05, 
consumed.food.per.person=400, 
intensity.of.agriculture=1,
x=1.3595
)

#Cambio 1
parameters2<-c(
growth.population.rate=0.0017, 
emigration.ratio=0.08, 
consumed.food.per.person=400, 
intensity.of.agriculture=0.8, 
x=1.3595)

#Cambio 2
parameters3<-c(
growth.population.rate=0.0017, 
emigration.ratio=0.1, 
consumed.food.per.person=400, 
intensity.of.agriculture=1.5,
x=1.3595)

InitialConditions <- c(
  population=100000, 
  fertility.of.agricultural.land=5000000,
  agricultural.land=8, 
  forest.land=5000
  )

times <- seq(-1000,2000,1)

# Definimos método de integración
 intg.method<-c("rk4")

#Especificamos modelo
MAYAS <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {

#auxiliary variables
food.produced <- agricultural.land*fertility.of.agricultural.land
demand.food <- population*consumed.food.per.person
food.gap <- demand.food - food.produced 
people.without.food <- food.gap/consumed.food.per.person

#flow variables
natural.population.increase <- growth.population.rate * population
emigration <- emigration.ratio*people.without.food
fertility.losses <- fertility.of.agricultural.land* min(2,(agricultural.land/forest.land)^x)/intensity.of.agriculture
deforestation <-min(food.gap/max(fertility.of.agricultural.land,1),forest.land/4)/intensity.of.agriculture
  
#state variables
dpopulation <- natural.population.increase - emigration
dfertility.of.agricultural.land <- (-1)*fertility.losses
dagricultural.land <- deforestation
dforest.land <- (-1)*deforestation

list(c(dpopulation, dfertility.of.agricultural.land, dagricultural.land, dforest.land))
    })
}

out1 <- ode(y = InitialConditions,
           times = times,
           func = MAYAS,
           parms = parameters1,
           method =intg.method)

out2 <- ode(y = InitialConditions,
           times = times,
           func = MAYAS,
           parms = parameters2,
           method =intg.method)
out3 <- ode(y = InitialConditions,
           times = times,
           func = MAYAS,
           parms = parameters3,
           method =intg.method)


datos1 <- data.frame(out1)
datos2 <- data.frame(out2)
datos3 <- data.frame(out3)

library(ggplot2)

#Gráfica de población 
 ggplot() +               
  geom_line(data = datos1, aes(time,population), 
             color = "dark green")+
  geom_line(data = datos2, aes(time,population), 
             color = "orange")+
    geom_line(data = datos3, aes(time,population), 
           color = "blue")+
  labs(x = "Time (years)", y = "Mayan population (persons)")+
  ggtitle("Mayan population") +  theme_minimal()+
  theme(plot.title = element_text(hjust = 0.5))

#Gráfica de fertiliy of agricultural land 
 ggplot() +               
  geom_line(data = datos1, aes(time,fertility.of.agricultural.land), 
             color = "darkgreen")+
   
  geom_line(data = datos2, aes(time,fertility.of.agricultural.land), 
             color = "orange")+
   
    geom_line(data = datos3, aes(time,fertility.of.agricultural.land), 
           color = "blue")+ 
   
  labs(x = "Time (years)", y = "Fertility of agricultural land (kg/km2)")+
  ggtitle("Behaviour of fertility of agricultural land (kg/km2)") +  theme_minimal()+
  theme(plot.title = element_text(hjust = 0.5))

 #Gráfica de  agricultural land 
 ggplot() +               
  geom_line(data = datos1, aes(time,agricultural.land), 
             color = "darkgreen")+
   
  geom_line(data = datos2, aes(time,agricultural.land), 
             color = "orange")+
   
    geom_line(data = datos3, aes(time,agricultural.land), 
           color = "blue")+ 
   
  labs(x = "Time (years)", y = " Agricultural land (km2)")+
  ggtitle("Behaviour of  agricultural land (km2)") +  theme_minimal()+
  theme(plot.title = element_text(hjust = 0.5))

 #Gráfica de  forest land 
 ggplot() +               
  geom_line(data = datos1, aes(time,forest.land), 
             color = "darkgreen")+
   
  geom_line(data = datos2, aes(time,forest.land), 
             color = "orange")+
   
    geom_line(data = datos3, aes(time,forest.land), 
           color = "blue")+ 
   
  labs(x = "Time (years)", y = " Forest land (km2)")+
  ggtitle("Behaviour of forest land (km2)") +  theme_minimal()+
  theme(plot.title = element_text(hjust = 0.5))

Descripción de la gráfica: El color verde de la gráfica representa el comportamiento de las variables de estado del sistema con un consume food per person 400 kg por persona al año, el color naranja representa el comportamiento con un consume food per person 200 kg por persona al año, mientras que la línea de color azul representa el comportamiento con 600 kg por persona al año, de consume food per person. Con base a las gráficas podemos decir que el sistema no es sensible a cambios marginales en los párametros.

5. Propón una política que evite el colapso de la civilización maya. Implementa esta política en el modelo.

Para evitar el colapso de la población se decició implementar la política pública denominada “Fertilization workshop”, la cual impacta directamente a la variable de estado fertily of agricultural land lo cual permite que la demand of food no crezca más que la food produced y por tanto la población maya no tenga que emigrar. Todo lo anterior se puede observar en el CLD presentado a conitnuación, así como el comportamiento de la población (en color morado) que no muestra un colpaso a diferencia del comportamiento representado en color azul, donde la population colapsa en el año 800 A.D.

#With political
parameters<-c(
growth.population.rate=0.0017,  
emigration.ratio=0.05, 
consumed.food.per.person=400, 
intensity.of.agriculture=1,
x=1.3595,
fertilization.per.km= 1 #km2
)

 InitialConditions <- c(
#Siempre debe serguir el mismo orden cuando se definen.
  population=100000, 
  fertility.of.agricultural.land=5000000,
  agricultural.land=8, 
  forest.land=5000
  )

times <- seq(-1000,2000,1)

# Definimos método de integración
 intg.method<-c("rk4")

#Especificamos modelo
MAYAS <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {

#auxiliary variables
food.produced <- agricultural.land*fertility.of.agricultural.land
demand.food <- population*consumed.food.per.person
food.gap <- demand.food - food.produced 
people.without.food <- food.gap/consumed.food.per.person

#flow variables
natural.population.increase <- growth.population.rate * population
emigration <- emigration.ratio*people.without.food
fertility.losses <- fertility.of.agricultural.land* min(2,(agricultural.land/forest.land)^x)/intensity.of.agriculture
deforestation <-min(food.gap/max(fertility.of.agricultural.land,1),forest.land/4)/intensity.of.agriculture
fertilizers <- population*fertilization.per.km
  
#state variables
dpopulation <- natural.population.increase - emigration
dfertility.of.agricultural.land <-  fertilizers- fertility.losses
dagricultural.land <- deforestation
dforest.land <- (-1)*deforestation

list(c(dpopulation, dfertility.of.agricultural.land, dagricultural.land, dforest.land))
    })
}

out.political <- ode(y = InitialConditions,
           times = times,
           func = MAYAS,
           parms = parameters,
           method =intg.method)

data.political <- data.frame(out.political)

#Gráfica de población 
 population.political <- ggplot (data.political, aes(x = time, y = population)) + 
   geom_line(color="purple")+
  xlab('Time (years)') + 
  ylab('Mayan Population (persons)') +
  ggtitle('Behaviuor population with Policy') + 
  theme_minimal()+
  theme(plot.title = element_text(hjust = 0.5))
 
 
 library(ggpubr)

comparation <- ggarrange(population.collapse, population.political,
                     ncol=2,nrow=1
                     )
comparation